Mean-field dynamics to negative absolute temperatures in the Bose-Hubbard model 
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We apply time-dependent Gutzwiller mean-field theory to provide a qualitative understanding for 
bosons in optical lattices that approach states corresponding to negative absolute temperatures. We 
perform the dynamical simulations to relate to the recent experiments by Braun et al, Ref. [J]. Time- 
of-flight images calculated from the two-dimensional numerical simulations reproduce characteristics 
of the experimental observations, in particular, the emergence of the four peaks at the corners of 
the Brillouin zone. 
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Physical systems with a perfect isolation and an upper 
bound on the energy are allowed to be in equilibrium 
states that are described by a negative absolute tem- 
perature, T < 0. These conditions are not too severe, 
yet such systems are only rarely observed. The main 
reason is that most systems are described by Hamiltoni- 
ans that have no upper bounds: free particles with mass 
M > have kinetic energies p 2 /2M, or using a more 
specific example, ultracold atoms are trapped by poten- 
tials V(r) rj V r 2 with Vq > 0. Negative T was first 
observed in nuclear spin systems [U, El] , which have very 
weak coupling to the environment, and naturally possess 
both lower and upper bounds on the energy due to the 
finite-dimensional Hilbert space of the spins. It is much 
harder to achieve negative absolute temperatures for sys- 
tems with motional degrees of freedom. Nevertheless, ul- 
tracold atoms in optical lattices hold the possibility to 
observe this novel state of matter [H, Q, since in these 
atomic clouds the system parameters can be controlled 
experimentally with great flexibility and accuracy @. 

In a recent experiment, Braun et al. [l[ indeed showed 
that attractively interacting bosonic atoms in an anti- 
trapping potential with Vq < can relax to a ther- 
modynamically stable state with a macroscopic occupa- 
tion (BEC) of the highest-energy single-particle states. 
This is also a demonstration of the inverse population 
of energy levels, a striking property of systems with 
T < 0, which can lead to counter-intuitive phenom- 
ena. Nevertheless, ultracold atomic experiments always 
need to start at T > and Vq > 0, and a certain 
protocol (see, e.g., Ref. @) is required to reach the 
Vq < side. Although numerical simulations for two- 
component fermionic clouds [f| were performed to have 
an understanding about the dynamics as the T < state 
is approached, bosons were not yet studied in a time- 
dependent simulation. In this work we address this ques- 
tion to complement the results for fermions and, more im- 
portantly, for comparison with the current experiments. 

Experimental setup and model. We are modeling 
the experiments discussed in Ref. [l|, where 39 K atoms 
are used to reach T < 0, based on the protocol outlined 
below. In 25 ms, a three-dimensional optical lattice at 
laser wavelength Al — 736.65nm is ramped up to a lat- 
tice depth s = 22E R , where E R = h 2 k\/{2M) is the 



recoil energy and kr, — 2tt/\t,. Using this deep lattice, 
in combination with a large positive scattering length 
a s = 309aBoiu- and a tight harmonic confinement, most 
of the cloud is in a Mott phase with one atom per lat- 
tice site. In the next step, in 2 ms, the magnetic field 
is ramped over a Feshbach resonance to reach a nega- 
tive scattering length, a s < 0, and simultaneously the 
red-detuned dipole traps in the horizontal directions are 
also ramped down. Thus the blue-detuned optical lat- 
tice beams provide an effective anti-trapping potential, 
Vq < 0. After waiting for 1 ms, the horizontal opti- 
cal lattice beams are ramped down to the final lattice 
depth Shor,/ = 6-Er in 2.5 ms. Due to technical difficul- 
ties regarding the compensation of gravity and to reach 
strong anti-trapping potentials, the vertical optical lat- 
tice is kept at a depth of s vcr = 22Er, and the tight 
vertical confinement is never reversed. We will use this 
detail in our advantage to greatly simplify the numerical 
simulations. 

With the assumption that the atoms are confined to 
the lowest Bloch band of the lattice 6], they can be de- 
scribed by the Bose-Hubbard model 7,Q, defined by 
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where (ij) denotes nearest neighbors, J is the nearest 
neighbor hopping rate, Vq characterizes the strength of 
the harmonic potential, [1q is a chemical potential to fix 
the number of atoms and U is the on-site interaction 
strength. Bosons are created and destroyed at a site j by 
bj and bp respectively, and hj = b^by 

Approximately 10 5 atoms were used in the three- 
dimensional setup of the experiment. We simplify our nu- 
merical calculations and consider only a two-dimensional 
horizontal layer of the atomic cloud. The simulations 
start at t = 20 ms of the experiment. At this time, the 
lattice depth is already quite deep, s sa 17. 6Er, and we 
neglect the vertical hopping between the layers. 

For simplicity, we calculate the hopping J and interac- 
tion strength U in harmonic approximation, where they 
depend on the lattice depths and the magnetic field B 
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U(s hor ,s VCI ,B)/E R = J — k L a s (B) s^s^J. (3) 

In the case of the interaction, the anisotropy of the Wan- 
nier functions in the horizontal and vertical directions 
has been taken into account. For 39 K atoms the scatter- 
ing length as a function of the magnetic field B near the 
Feshbach resonance used in the experiment is approxi- 
mately @ 



a s (B)=a hg [1-A/(B-B les )\, 



(4) 



where ab g = — 33asohr, A = — 52G, and B ICS = 402. 4G. 

Method. To simulate the real-time dynamics of 
bosons, we apply time-dependent Gutzwiller mean- field 
theory. The bosonic Gutzwiller ansatz (GA) was first 
applied to approximate the ground state of the homoge- 
neous version ofEq. (TJ [T3, El- The time-dependent G A 
for inhomogeneous systems wa s appli ed for various non- 
equilibrium scenarios recently |12h16| . In this approach, 
the time evolution of the m-boson occupation amplitude 
at site j of the lattice and time t, denoted by f m {j, t), is 
given by the following differential equation: 
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-J(t) </>*(j,t)Vrri + lf m +i(j,t) 
-J{t)4> (j,t) V^/m-i(j,*). (5) 



The mean-field decoupling is 
the hopping term of Eq. (fTJ), 
.sJV^+T /* (j + -5, t)f m+1 (j + S, 1)] 



applied in 

W,t) = 

E«,m[>/m+T/*(j+<y,*)/flH-i(j +/,*)]> wh ere S 
runs over the nearest neighbors of j. The amplitudes are 
normalized, 1 = J2 m \fm(j,t)\ 2 , Vt,j. 

The explicit time-dependence of the microscopic pa- 
rameters U(t) and J(t) are computed using Eqns. @, ^ 
and ([4]) from the time-dependent optical lattice depths 
Shor(i) and s vcr t(t) and the magnetic field B(t), while 
Vo(t) is determined from the horizontal frequency using 
Vo(t) ~ Whor(*)- More details can be found in Ref. 

At to = 20 ms, the initial condition f m (j,to) is de- 
termined in two steps. At each site the amplitudes are 
calculated first in GA in the local density approxima- 
tion, followed by a few self-consistency sweeps until the 
stationary solution of Eq. (JSJ) is found. We used square 
lattices between 80 x 80 and 160 x 160 lattice sites and 
an initial atom number of iVt t(£o) ~ 1920. The maximal 
allowed occupation number was usually m c = 6. 

To numerically integrate Eq. ([5]), we use a split-step 
method. At each time step t — >• t + St, we first perform 
unitary rotations corresponding to the diagonal terms in 
Eq. ([5]). Second, we update the mean-field <j). Third, we 



perform an Euler step using the remaining terms. As the 
last step, we normalize the amplitudes. 

Next, we analyze the results of the simulations and 
compare them to the experiments. The naming conven- 
tions of the various protocols and the corresponding final 
experimental parameters are summarized in Table U 



Protocol name 


Uf ~ a sJ 


Voj ~ (uj 1ioi j/2tv) 2 


a) 


— 37 dBohr 


- 43 2 sec~ 2 


c) 


— 37 aBohr 


+ 35 2 sec -2 


d) 


— 37 aBohr 


+ 42 2 sec~ 2 


e) 


— 37 aBohr 


+ 60 2 sec -2 


0) 


+ 33 aBohr 


+ 44 2 sec- 2 



TABLE I: Naming convention for the different protocols. The 
first parts of the protocols (t < 25 ms) are identical, the 
difference is in the final values of the interaction strengths 
and the external harmonic potentials. See also Figs. 1 1 121 

Results. We calculate various macroscopic quantities 
using the simulation, the total atom number 



AW*) = n i^' 



with rij(t) = J2 m m \fm(j, t)\ 2 , the cloud radius 
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R W = ^rfoityNUt), 

the condensate fraction 

N (t) = El^W>l 2 ' 



(7) 



(8) 



with (bj(t)) = E m V^/;(i.*)/m+i(j^)» and the 
nearest neighbor superfluid correlations 

C(t) = ^{btitM)) ^ (9) 

<ij> <ij> 

where the last equation follows from the GA. Note that 
C is real and related to the kinetic energy, given by 



Kit) = -J(t)C(t) 



(10) 



We plot these quantities along with microscopic param- 
eters for different protocols in Figs. [T] and [2] Although 
iVtot is conserved by Eq. ((SJ), there is a weak time depen- 
dence due to errors of the numerical integration. 

We see that when the signs of the final values of the 
interactions, Uf, and the harmonic potential, Vqj, are 
the same, the system approaches a macroscopically sta- 
tionary state. This is consistent both with the theoretical 
expectations that the cloud can equilibrate when the en- 
ergy spectrum is bounded from at least one side, and 
with the experimental demonstration of the thermody- 
namic stability in these cases [l[. There is, however, a 
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FIG. 1: (color online) Microscopic parameters (upper panels) 
and the corresponding macroscopic quantities (lower panels) 
from time-dependent Gutzwiller simulations for the protocols 
a) (left panels) and 0) (right panels). The system approaches 
a macroscopically stationary state even for attractive final 
interactions Uf < and anti-trapping potentials Voj < 0. 



striking difference between the repulsively and the at- 
tractively interacting cases. In the more usual, repul- 
sive case, with Uf > and Voj > 0, the average super- 
fluid correlations on nearest neighbor sites are positive, 
C > 0, and the system has a negative total kinetic en- 
ergy. In contrast, in the attractive case, with Uf < 
and Voj < , the superfluid order parameter on near- 
est neighbors shows, on average, anticorrelation, C < 0, 
which directly translates to a positive total kinetic en- 
ergy, cf. Eq. (fit)]) . With the symmetric single-particle 
band and weak interactions, K > can only occur with 
an inverse population of the single-particle energy levels. 
From a different perspective, the SF order parameter is 
alternating on the two sublattices, (bj) w (— iy (bo) in the 
case a). Thus the spatial Fourier transform (related also 
to time-of-flight images, see below) is enhanced around 
momenta Q = (±fcL,±fci). We note that the ground 
state of Eq. ^ with parameters sgn([7 f) J, \ Uf\ and \ Voj\ 
and the same N to t as in the protocols a) or 0) have su- 
perfluid fractions /N to t ~ 0.97 and cloud radiuses 
i? (0) /-R(20ms) m 1.1 in GA. These values suggest that 
the stationary states are "heated" in comparison to the 
corresponding ground states. 

Quite remarkably, in addition to the protocols which 
lead to thermodynamically stable clouds in experiments, 
the Gutzwiller approach also seems to capture qualita- 
tively correctly the cases with attractive final interac- 
tions, Uf < and trapping potentials, Voj > 0. In these 
cases, the Hamiltonian is not bounded, and equilibrium 
is not possible at any finite 1/T ^ 0. The dynamics is 
governed by two microscopic processes: the cloud is al- 
lowed to expand (see Fig. [2]) , provided that the increase 
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FIG. 2: (color online) Macroscopic quantities from time- 
dependent Gutzwiller simulations for the protocols c) (red 
solid line), d) (blue dashed line) and e) (green dashed-dotted 
line). The system expands and does not relax to a stationary 
state. The SF coherence decays on short time scales. The 
vertical dashed lines correspond to t = 30.5 ms. Data for 
protocol a) is shown for comparison (black dotted line). 



in potential energy is covered by a decrease in interac- 
tion energy (not shown), i.e., the atoms try to cluster. 
However, the expansion also leads to the dilution of the 
cloud, and thus the probability of clustering decreases. 
Therefore the expansion gets slower and slower. 

Finally, we focus on time-of-flight (TOF) images. Fol- 
lowing Ref . @ , we define the (two-dimensional) TOF in- 
tensity as 



iroF(k) ~ Hk)| 2 S(k) 
with the single-particle correlation function 



(11) 



,ik(ri-rj) 



Hbtb )^N tot -No + \b(K)\\ (12) 



where 6(k) = J2 r 



-?'kr, 



(bj). The Fourier transform of 



the Wannier function in harmonic approximation gives 



Kk)|< 



-k 2 /(v^fe|) 



We plot TOF images in 



Gutzwiller approximation for the different protocols and 
different times in Fig. [3J to be compared with the similar 
images in Ref. [l[. Based on the TOF images, we also 
define visibilities V = (Ia — Ib)/{Ia + Ib), with Ia,b 
the intensities integrated in a small area around vectors 
Q = [Ul, and Q = (v2fcc, 0), respectively. These are 
shown in Fig. 21 and we see a fast decay to negative visi- 
bilities, similar to the experiments. However, we were not 
able to extract lifetimes, mainly due to the fluctuations 
in the data which can come from finite size effects. 

Conclusions. We used time-dependent Gutzwiller 
mean-field theory to describe the dynamics of bosons 
which arc subject to various protocols with the aim of 
achieving T < with ultracold atoms in optical lattices. 
We found that these numerical results qualitatively agree 
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FIG. 3: (color online) Two-dimensional time-of-flight images from the Gutzwiller calculation for the different protocols at 
different times. The intensities Ttof (k) have been normalized by the total particle number. 
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FIG. 4: (color online) Visibilities (defined similarly as in 
Ref. 01) as a function of time t. For protocol a), we see no 
decay. For the protocols c), d) and e), the visibilities decay 
rapidly to negative final values, similar to the experiments. 

with the experimental findings. Most importantly, after 
reversing the signs of both the interactions between the 
atoms and the external harmonic potential, the system 



relaxes to a macroscopically stationary state where the 
BEC shows up at finite momenta. This novel BEC pat- 
tern can be directly seen in the TOF images. The simple 
time-dependent mean-field theory also seems to capture 
the experimental observation that the combination of an 
attractive final interaction and a trapping final potential 
leads to a rapid decay in TOF images. 

Despite that we found qualitative similarities for some 
quantities between the numerics and the experiments, 
there are limitations which prohibit further, quantita- 
tive agreement. The experiment operates on a three- 
dimensional cloud, which is too large for the current nu- 
merics. The evaluation of microscopic parameters can 
also be improved by using the integrals of the numerically 
exact Wannier functions. Additionally, thermodynamic 
quantities, like entropy or temperature, cannot be intro- 
duced in this GA in a natural way. Addressing these and 
further restrictions are beyond the scope of the present 
work and left for future studies. 
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